Section 1 Conditions environnementales régionales

Est-ce que les régions varient dans leurs conditions environnementales?

#Beaufort Sea: Mackenzie Shelf, Amundsen Gulf Maud
#Kitikmeot: Coronation Maud, Larsen Sound, Peel Sound
#Baffin Bay: North Water, Lancaster Sound, NW Baffin Bay
dfish$province <- with(dfish, ifelse(region %in% c("Mackenzie Shelf", "Amundsen Gulf Mouth"), "Beaufort Sea",
                                    ifelse(region %in% c("Coronation Maud", "Larsen Sound - Victoria Strait", "Peel Sound"), "Kitikmeot",
                                           ifelse(region == "NEG", "Greenland Sea",
                                                  "NW Baffin Bay"))))

dfish$region_year <- as.factor(paste(dfish$region, dfish$year))
dfish$region_year <- ordered(dfish$region_year, levels = c("Mackenzie Shelf 2009","Mackenzie Shelf 2010","Mackenzie Shelf 2014","Mackenzie Shelf 2015",
                                                           "Amundsen Gulf Mouth 2015", 
                                                           "Coronation Maud 2011" ,"Coronation Maud 2016","Coronation Maud 2017","Coronation Maud 2018",
                                                           "Larsen Sound - Victoria Strait 2010", "Larsen Sound - Victoria Strait 2016","Larsen Sound - Victoria Strait 2018",
                                                           "Peel Sound 2010","Peel Sound 2011","Peel Sound 2014","Peel Sound 2015","Peel Sound 2016",
                                                           "Lancaster Sound 2010","Lancaster Sound 2011","Lancaster Sound 2015" ,"Lancaster Sound 2016","Lancaster Sound 2017",
                                                           "North Water 2014","North Water 2016", "North Water 2018",
                                                           "West Baffin Bay 2015","West Baffin Bay 2016",
                                                           "NEG 2017"
                                                           
                                                           ))
station_data <- aggregate(data = dfish, unique_fish_id ~ region_year + station + surf_temp_degC + surf_sal_kgm3 + open_water_day + NASC_zoo + indstandard + prof_mel, length)
env_pca <- station_data%>% select(surf_temp_degC, surf_sal_kgm3, open_water_day, region_year, prof_mel) %>% PCA(quali.sup = c(4), graph = FALSE)

station_data$pc1 <- env_pca$ind$coord[, 1]
station_data$pc2 <- env_pca$ind$coord[, 2]

pca.vars <- env_pca$var$coord %>% data.frame
pca.vars$vars <- rownames(pca.vars)

posArrowsX <- c(4, 4, 4,4)
posArrowsY<- c(4,4,4,4)
colours <- c(pals::stepped(n=24)[13:16], pals::polychrome(n=4)[4], pals::stepped(n=24)[5:11], pals::stepped2(n=20)[13], pals::stepped(n=24)[c(1:4,17:20)],pals::stepped3(n=16)[16],pals::stepped(n=24)[21:23], pals::stepped3(n=20)[c(5:6)], "red")
names(colours) <- levels(station_data$region_year)
gg_PCA_env <- ggplot() +  
  geom_point(data = station_data, aes(x = pc1, y = pc2, colour = region_year), size = 2.5)+
  #scale_colour_manual(values = colours)+
  coord_equal()+
  theme_classic()+
  theme(panel.background = element_rect(fill = "transparent"), # bg of the panel
        plot.background = element_rect(fill = "transparent", color = NA),
        legend.background = element_rect(fill=alpha('white', 0.4)))+
  labs(x= paste("PC1 (", round(env_pca$eig["comp 1", "percentage of variance"], digits = 2), "%)"),
       y= paste("PC2 (", round(env_pca$eig["comp 2", "percentage of variance"], digits = 2), "%)")) +
  geom_hline(yintercept = 0, lty = 2) +
  geom_vline(xintercept = 0, lty = 2) +
  geom_segment(data = pca.vars, aes(x = 0, xend = Dim.1*3, y = 0, yend = Dim.2*3),
               arrow = arrow(length = unit(0.025, "npc"), type = "open"), 
               lwd = 0.5) + 
  geom_text(data = pca.vars, 
            aes(x = Dim.1*posArrowsX, y =  Dim.2*posArrowsY, label = vars), 
            check_overlap = FALSE)+
  scale_color_manual(values = colours)
gg_PCA_env  

envzoo_pca <- station_data%>% select(surf_temp_degC, surf_sal_kgm3, open_water_day, NASC_zoo, region_year, prof_mel) %>% PCA(quali.sup = c(5), graph = FALSE)

station_data$pc1_zoo <- envzoo_pca$ind$coord[, 1]
station_data$pc2_zoo <- envzoo_pca$ind$coord[, 2]

pca.vars <- envzoo_pca$var$coord %>% data.frame
pca.vars$vars <- rownames(pca.vars)

posArrowsX <- c(4, 4, 4,4,4)
posArrowsY<- c(4,4,4,4,4)
gg_PCA_envzoo <- ggplot() +  
  geom_point(data = station_data, aes(x = pc1_zoo, y = pc2_zoo, colour = region_year), size = 2.5)+
  #scale_colour_manual(values = colours)+
  coord_equal()+
  theme_classic()+
  theme(panel.background = element_rect(fill = "transparent"), # bg of the panel
        plot.background = element_rect(fill = "transparent", color = NA),
        legend.background = element_rect(fill=alpha('white', 0.4)))+
  labs(x= paste("PC1 (", round(envzoo_pca$eig["comp 1", "percentage of variance"], digits = 2), "%)"),
       y= paste("PC2 (", round(envzoo_pca$eig["comp 2", "percentage of variance"], digits = 2), "%)")) +
  geom_hline(yintercept = 0, lty = 2) +
  geom_vline(xintercept = 0, lty = 2) +
  geom_segment(data = pca.vars, aes(x = 0, xend = Dim.1*3, y = 0, yend = Dim.2*3),
               arrow = arrow(length = unit(0.025, "npc"), type = "open"), 
               lwd = 0.5) + 
  geom_text(data = pca.vars, 
            aes(x = Dim.1*posArrowsX, y =  Dim.2*posArrowsY, label = vars), 
            check_overlap = FALSE)+
  scale_color_manual(values = colours)
gg_PCA_envzoo

envfish_pca <- station_data%>% select(surf_temp_degC, surf_sal_kgm3, open_water_day, indstandard, region_year, prof_mel) %>% PCA(quali.sup = c(5), graph = FALSE)

station_data$pc1_fish <- envfish_pca$ind$coord[, 1]
station_data$pc2_fish <- envfish_pca$ind$coord[, 2]

pca.vars <- envfish_pca$var$coord %>% data.frame
pca.vars$vars <- rownames(pca.vars)

posArrowsX <- c(4, 4, 4,4,4)
posArrowsY<- c(4,4,4,4,4)
gg_PCA_envfish <- ggplot() +  
  geom_point(data = station_data, aes(x = pc1_fish, y = pc2_fish, colour = region_year), size = 2.5)+
  #scale_colour_manual(values = colours)+
  coord_equal()+
  theme_classic()+
  theme(panel.background = element_rect(fill = "transparent"), # bg of the panel
        plot.background = element_rect(fill = "transparent", color = NA),
        legend.background = element_rect(fill=alpha('white', 0.4)))+
  labs(x= paste("PC1 (", round(envfish_pca$eig["comp 1", "percentage of variance"], digits = 2), "%)"),
       y= paste("PC2 (", round(envfish_pca$eig["comp 2", "percentage of variance"], digits = 2), "%)")) +
  geom_hline(yintercept = 0, lty = 2) +
  geom_vline(xintercept = 0, lty = 2) +
  geom_segment(data = pca.vars, aes(x = 0, xend = Dim.1*3, y = 0, yend = Dim.2*3),
               arrow = arrow(length = unit(0.025, "npc"), type = "open"), 
               lwd = 0.5) + 
  geom_text(data = pca.vars, 
            aes(x = Dim.1*posArrowsX, y =  Dim.2*posArrowsY, label = vars), 
            check_overlap = FALSE)+
  scale_color_manual(values = colours)
gg_PCA_envfish

envbio_pca <- station_data%>% select(surf_temp_degC, surf_sal_kgm3, open_water_day, NASC_zoo,indstandard, region_year, prof_mel) %>% PCA(quali.sup = c(6), graph = FALSE)

station_data$pc1_bio <- envbio_pca$ind$coord[, 1]
station_data$pc2_bio <- envbio_pca$ind$coord[, 2]

pca.vars <- envbio_pca$var$coord %>% data.frame
pca.vars$vars <- rownames(pca.vars)

posArrowsX <- c(4, 4, 4,4,4,4)
posArrowsY<- c(4,4,4,4,4,4)
gg_PCA_bio <- ggplot() +  
  geom_point(data = station_data, aes(x = pc1_bio, y = pc2_bio, colour = region_year), size = 2.5)+
  #scale_colour_manual(values = colours)+
  coord_equal()+
  theme_classic()+
  theme(panel.background = element_rect(fill = "transparent"), # bg of the panel
        plot.background = element_rect(fill = "transparent", color = NA),
        legend.background = element_rect(fill=alpha('white', 0.4)))+
  labs(x= paste("PC1 (", round(envbio_pca$eig["comp 1", "percentage of variance"], digits = 2), "%)"),
       y= paste("PC2 (", round(envbio_pca$eig["comp 2", "percentage of variance"], digits = 2), "%)")) +
  geom_hline(yintercept = 0, lty = 2) +
  geom_vline(xintercept = 0, lty = 2) +
  geom_segment(data = pca.vars, aes(x = 0, xend = Dim.1*3, y = 0, yend = Dim.2*3),
               arrow = arrow(length = unit(0.025, "npc"), type = "open"), 
               lwd = 0.5) + 
  geom_text(data = pca.vars, 
            aes(x = Dim.1*posArrowsX, y =  Dim.2*posArrowsY, label = vars), 
            check_overlap = FALSE)+
  scale_color_manual(values = colours)
gg_PCA_bio

hist(aggregate(open_water_day ~ year + region, data = dfish, FUN = mean)$open_water_day)

plot(open_water_day ~ region, data = dfish, las = 2)

plot(open_water_day ~ as.factor(province), data = dfish)

aggregate(open_water_day ~ year + region +sampling_date + bu_date, data = dfish, FUN = mean)
##    year                         region sampling_date    bu_date open_water_day
## 1  2009                Mackenzie Shelf    2009-08-05 2009-06-22             44
## 2  2009                Mackenzie Shelf    2009-08-15 2009-06-22             54
## 3  2010                Mackenzie Shelf    2010-08-16 2010-06-07             70
## 4  2010                Mackenzie Shelf    2010-08-19 2010-06-07             73
## 5  2010                Mackenzie Shelf    2010-08-21 2010-06-07             75
## 6  2010                Mackenzie Shelf    2010-08-22 2010-06-07             76
## 7  2010                Mackenzie Shelf    2010-08-23 2010-06-07             77
## 8  2010                Mackenzie Shelf    2010-08-24 2010-06-07             78
## 9  2010                Lancaster Sound    2010-08-07 2010-07-05             33
## 10 2010                Lancaster Sound    2010-08-08 2010-07-05             34
## 11 2010 Larsen Sound - Victoria Strait    2010-08-10 2010-08-09              1
## 12 2010                     Peel Sound    2010-08-09 2010-08-16             -7
## 13 2011                Lancaster Sound    2011-08-04 2011-05-23             73
## 14 2011                Coronation Maud    2011-08-10 2011-07-18             23
## 15 2011                     Peel Sound    2011-08-08 2011-07-25             14
## 16 2014                    North Water    2014-08-01 2014-05-19             74
## 17 2014                    North Water    2014-08-02 2014-05-19             75
## 18 2014                    North Water    2014-08-06 2014-05-19             79
## 19 2014                Mackenzie Shelf    2014-08-23 2014-06-23             61
## 20 2014                Mackenzie Shelf    2014-08-24 2014-06-23             62
## 21 2014                     Peel Sound    2014-08-10 2014-09-01            -22
## 22 2014                     Peel Sound    2014-08-11 2014-09-01            -21
## 23 2015            Amundsen Gulf Mouth    2015-08-24 2015-05-11            105
## 24 2015                Lancaster Sound    2015-08-10 2015-05-18             84
## 25 2015                Lancaster Sound    2015-08-11 2015-05-18             85
## 26 2015                Lancaster Sound    2015-08-13 2015-05-18             87
## 27 2015                Mackenzie Shelf    2015-08-27 2015-05-25             94
## 28 2015                Mackenzie Shelf    2015-08-30 2015-05-25             97
## 29 2015                West Baffin Bay    2015-08-06 2015-07-06             31
## 30 2015                West Baffin Bay    2015-08-07 2015-07-06             32
## 31 2015                     Peel Sound    2015-08-15 2015-08-10              5
## 32 2016                    North Water    2016-08-04 2016-05-16             80
## 33 2016                    North Water    2016-08-06 2016-05-16             82
## 34 2016                    North Water    2016-08-07 2016-05-16             83
## 35 2016                    North Water    2016-08-08 2016-05-16             84
## 36 2016                    North Water    2016-08-09 2016-05-16             85
## 37 2016                    North Water    2016-08-10 2016-05-16             86
## 38 2016                Lancaster Sound    2016-08-17 2016-05-23             86
## 39 2016                Lancaster Sound    2016-08-18 2016-05-23             87
## 40 2016                West Baffin Bay    2016-08-01 2016-07-04             28
## 41 2016                Coronation Maud    2016-08-21 2016-07-11             41
## 42 2016                Coronation Maud    2016-08-22 2016-07-11             42
## 43 2016                Coronation Maud    2016-08-23 2016-07-11             43
## 44 2016 Larsen Sound - Victoria Strait    2016-08-20 2016-07-25             26
## 45 2016                     Peel Sound    2016-08-19 2016-08-01             18
## 46 2017                Lancaster Sound    2017-08-03 2017-05-22             73
## 47 2017                Lancaster Sound    2017-08-04 2017-05-22             74
## 48 2017                Lancaster Sound    2017-08-05 2017-05-22             75
## 49 2017                            NEG    2017-08-25 2017-06-19             67
## 50 2017                            NEG    2017-08-29 2017-06-19             71
## 51 2017                            NEG    2017-08-31 2017-06-19             73
## 52 2017                            NEG    2017-09-02 2017-06-19             75
## 53 2017                            NEG    2017-09-04 2017-06-19             77
## 54 2017                            NEG    2017-09-06 2017-06-19             79
## 55 2017                Coronation Maud    2017-08-08 2017-07-03             36
## 56 2017                Coronation Maud    2017-08-09 2017-07-03             37
## 57 2018                    North Water    2018-08-28 2018-05-21             99
## 58 2018                    North Water    2018-08-29 2018-05-21            100
## 59 2018                Coronation Maud    2018-08-21 2018-07-16             36
## 60 2018                Coronation Maud    2018-08-22 2018-07-16             37
## 61 2018 Larsen Sound - Victoria Strait    2018-08-19 2018-08-13              6
lattice::xyplot(open_water_day ~ as.factor(year)|region, data = dfish)

hist(aggregate(surf_temp_degC ~ year + station, data = dfish, FUN = mean)$surf_temp_degC)

plot(surf_temp_degC ~ region, data = dfish, las = 2) #Beaucoup de variation Mackenzie Self

aggregate(surf_temp_degC ~ year + region +sampling_date + station + NASC_zoo + indstandard + surf_sal_kgm3 + open_water_day + bu_date, data = dfish, FUN = mean) %>% filter(region == "Mackenzie Shelf") 
##    year          region sampling_date station  NASC_zoo indstandard
## 1  2009 Mackenzie Shelf    2009-08-05     260 12.175651 0.266204537
## 2  2009 Mackenzie Shelf    2009-08-15     345 12.175651 0.010588621
## 3  2010 Mackenzie Shelf    2010-08-16      14 34.703792 0.011038416
## 4  2010 Mackenzie Shelf    2010-08-16      10 34.703792 0.022063208
## 5  2010 Mackenzie Shelf    2010-08-16      12 34.703792 0.013204152
## 6  2010 Mackenzie Shelf    2010-08-16      15 34.703792 0.005333019
## 7  2010 Mackenzie Shelf    2010-08-19       5 34.703792 0.007695008
## 8  2010 Mackenzie Shelf    2010-08-21       1 34.703792 0.038703415
## 9  2010 Mackenzie Shelf    2010-08-21       2 34.703792 0.016440352
## 10 2010 Mackenzie Shelf    2010-08-22       6 34.703792 0.006108735
## 11 2010 Mackenzie Shelf    2010-08-22       8 34.703792 0.007410152
## 12 2010 Mackenzie Shelf    2010-08-23      16 34.703792 0.012712653
## 13 2010 Mackenzie Shelf    2010-08-24      13 34.703792 0.010047687
## 14 2014 Mackenzie Shelf    2014-08-23     434  4.490901 0.005448571
## 15 2014 Mackenzie Shelf    2014-08-24     421  4.490901 0.001690748
## 16 2015 Mackenzie Shelf    2015-08-27     435  3.203247 0.002783745
## 17 2015 Mackenzie Shelf    2015-08-30     421  3.203247 0.000604846
##    surf_sal_kgm3 open_water_day    bu_date surf_temp_degC
## 1       29.46308             44 2009-06-22      4.1338462
## 2       27.75806             54 2009-06-22      2.1182560
## 3       26.13440             70 2010-06-07      9.2080000
## 4       27.16917             70 2010-06-07      8.6690000
## 5       27.94775             70 2010-06-07      8.0842500
## 6       28.05033             70 2010-06-07      8.0592222
## 7       25.85400             73 2010-06-07      9.7882000
## 8       25.92320             75 2010-06-07     10.1243000
## 9       27.73389             75 2010-06-07      8.0745789
## 10      25.05264             76 2010-06-07      9.9577374
## 11      27.83875             76 2010-06-07      9.2170000
## 12      28.38987             77 2010-06-07      6.9733871
## 13      26.32015             78 2010-06-07      9.8555385
## 14      30.72773             61 2014-06-23      5.2609500
## 15      26.73137             62 2014-06-23      0.1061250
## 16      28.34078             94 2015-05-25      3.9685500
## 17      27.09729             97 2015-05-25     -0.3764118
plot(surf_temp_degC ~ as.factor(province), data = dfish)

lattice::xyplot(surf_temp_degC ~ as.factor(year)|region, data = dfish)

hist(aggregate(NASC_zoo ~ year+region, data = dfish, FUN = mean)$NASC_zoo)

plot(NASC_zoo ~ region, data = dfish, las = 2) #Beaucoup variation MS en 2010

plot(NASC_zoo ~ as.factor(province), data = dfish) #presque pas de zoo en mer du GL, et MS 2010  biaise vraiment la mer de Beaufort

plot(indstandard ~ region, data = dfish, las = 2) #Beaucoup variation LS -> c'est vraiment en 2010 qu'on a un pic

plot(log10(indstandard) ~ region, data = dfish, las = 2)

aggregate(indstandard ~ year + region +sampling_date + station, data = dfish, FUN = mean) %>% filter(region %in% c("Lancaster Sound")) 
##    year          region sampling_date   station indstandard
## 1  2010 Lancaster Sound    2010-08-07       301  0.90395209
## 2  2017 Lancaster Sound    2017-08-03       301  0.03870584
## 3  2017 Lancaster Sound    2017-08-04       304  0.01889668
## 4  2010 Lancaster Sound    2010-08-08       305  0.17822000
## 5  2016 Lancaster Sound    2016-08-18       305  0.02550317
## 6  2017 Lancaster Sound    2017-08-05       305  0.23590885
## 7  2011 Lancaster Sound    2011-08-04       332  0.09332530
## 8  2016 Lancaster Sound    2016-08-17 Allen Bay  0.03430469
## 9  2015 Lancaster Sound    2015-08-10     CAA-1  0.02836053
## 10 2015 Lancaster Sound    2015-08-10     CAA-2  0.11179698
## 11 2015 Lancaster Sound    2015-08-11     CAA-3  0.07652033
## 12 2015 Lancaster Sound    2015-08-13     CAA-4  0.10226513
## 13 2015 Lancaster Sound    2015-08-13     CAA-5  0.22170073
plot(surf_sal_kgm3 ~ region, data = dfish, las = 2) #Beaucoup variation CM & PS, oz on a d'ailleurs des salinités assez faibles

lattice::xyplot(surf_sal_kgm3 ~ as.factor(year)|region, data = dfish)

aggregate(surf_sal_kgm3 ~ year + region +sampling_date + station, data = dfish, FUN = mean) %>% filter(region %in% c("Coronation Maud", "Peel Sound")) #2011 PS, 2017 et 2018 CM (QMG1-2-3-M) -> semble dépendre où était située la station QMG -> plus loin vers l'ouest = plus salé
##    year          region sampling_date station surf_sal_kgm3
## 1  2014      Peel Sound    2014-08-10     309      30.46635
## 2  2010      Peel Sound    2010-08-09     310      28.32820
## 3  2011      Peel Sound    2011-08-08     310      19.93550
## 4  2014      Peel Sound    2014-08-11     310      30.28137
## 5  2016      Peel Sound    2016-08-19    310E      24.54000
## 6  2011 Coronation Maud    2011-08-10     314      25.85850
## 7  2016 Coronation Maud    2016-08-23     314      25.29824
## 8  2016 Coronation Maud    2016-08-23     316      27.13650
## 9  2015      Peel Sound    2015-08-15   CAA-7      29.13917
## 10 2016 Coronation Maud    2016-08-22   QMG-1      24.17040
## 11 2017 Coronation Maud    2017-08-09   QMG-1      21.01343
## 12 2018 Coronation Maud    2018-08-21   QMG-1      21.39480
## 13 2016 Coronation Maud    2016-08-21   QMG-2      20.82000
## 14 2017 Coronation Maud    2017-08-09   QMG-2      18.47175
## 15 2018 Coronation Maud    2018-08-21   QMG-2      23.47429
## 16 2016 Coronation Maud    2016-08-22   QMG-3      22.72938
## 17 2017 Coronation Maud    2017-08-09   QMG-3      27.54350
## 18 2016 Coronation Maud    2016-08-23   QMG-4      24.97867
## 19 2017 Coronation Maud    2017-08-09   QMG-4      27.97933
## 20 2018 Coronation Maud    2018-08-22   QMG-4      27.73485
## 21 2016 Coronation Maud    2016-08-21   QMG-M      19.80858
## 22 2017 Coronation Maud    2017-08-08   QMG-M      24.22665
## 23 2018 Coronation Maud    2018-08-22   QMG-M      25.57633
plot(prof_mel ~ region, data = dfish, las = 2) #Beaucoup variation CM & PS, oz on a d'ailleurs des salinités assez faibles

range(dfish$prof_mel)
## [1]   2.5 145.0
lattice::xyplot(prof_mel ~ as.factor(year)|region, data = dfish)

Section 2: Facteurs reliés aux poissons eux-mêmes

Quel est l’éventail de taille de mes poissons et quel impact la taille semble-t’elle avoir?

La taille des poissons varie entre les régions, avec des larves particulièrement grosses dans l’Amundsen Gulf Mouth et Mackenzie Shelf ainsi qu’un énorme individu (61mm) dans Peel Sound. Les larves les plus petites sont en général celle de West Baffin Bay, Victoria Strait (avec un accroissement) et Peel Sound.

Les poissons échantillonnés dans une zone de débâcle hâtive sont généralement un peu plus gros, bien qu’il y ait encore des petits dans le pool. Ça reste que la majorité des poissons sont très petits lorsqu’ils sont collectés à une période concordant à ou près de la débâcle.

La quantité de carbone ingéré semble augmenter avec la taille. Par contre, comme on a beaucoup moins de grosses larves dans le pool, le signal est quand même dur à détecter (on a un gros biais pour les petites larves).

Quand on regarde le succès alimentaire en fonction de la taille, le nuage de point est assez confus. Généralement, le succès alimentaire brut est faible, peu importe la taille. Quand on transforme avec log10, on voit qu’il y a une légère augmentation en fonction de la taille. La relation linéaire entre ces deux variables est significative. Par contre, la variabilité semble à son plus haut lorsque le poisson a une très petite taille et une grande taille. Donc, ça indiquerait que parmi les très petits, on a dû en attraper qui seraient morts de famine bientôt OU qui avaient encore des réserves du sac vitellin et les gros seraient peut-être plus résistants à une famine variable, d’où la variabilité. Entre 12-13 et 25mm, le succès alimentaire augmente avec moins de variabilité. Donc, bien qu’on corrige FS en divisant par la masse, la taille demeure un facteur important de capture de proie.

#taille des poissons
plot(est_standard_length ~ region, data = dfish, las = 2)

tapply(dfish$est_standard_length, INDEX = list(dfish$region, dfish$year), FUN = mean)
##                                   2009     2010     2011     2014     2015
## Amundsen Gulf Mouth                 NA       NA       NA       NA 32.85720
## Coronation Maud                     NA       NA 15.59700       NA       NA
## Lancaster Sound                     NA 13.23085 13.62750       NA 16.59357
## Larsen Sound - Victoria Strait      NA 11.30658       NA       NA       NA
## Mackenzie Shelf                12.6857 27.09488       NA 35.40000 31.54290
## NEG                                 NA       NA       NA       NA       NA
## North Water                         NA       NA       NA 17.73251       NA
## Peel Sound                          NA 10.81823 17.42556 12.13674 11.36489
## West Baffin Bay                     NA       NA       NA       NA 14.48783
##                                    2016     2017     2018
## Amundsen Gulf Mouth                  NA       NA       NA
## Coronation Maud                22.55682 16.26570 14.98442
## Lancaster Sound                17.96064 15.33450       NA
## Larsen Sound - Victoria Strait 14.74777       NA 11.03090
## Mackenzie Shelf                      NA       NA       NA
## NEG                                  NA 19.12333       NA
## North Water                    20.41569       NA 21.30636
## Peel Sound                     12.85020       NA       NA
## West Baffin Bay                11.19310       NA       NA
tapply(dfish$est_standard_length, INDEX = list(dfish$region, dfish$year), FUN = median)
##                                  2009   2010  2011   2014   2015   2016    2017
## Amundsen Gulf Mouth                NA     NA    NA     NA 35.000     NA      NA
## Coronation Maud                    NA     NA 16.07     NA     NA 18.906 16.1720
## Lancaster Sound                    NA 12.713 13.14     NA 15.857 18.438 14.6875
## Larsen Sound - Victoria Strait     NA 11.431    NA     NA     NA 14.063      NA
## Mackenzie Shelf                12.307 27.259    NA 38.000 33.500     NA      NA
## NEG                                NA     NA    NA     NA     NA     NA 18.6500
## North Water                        NA     NA    NA 18.924     NA 21.997      NA
## Peel Sound                         NA 10.790 12.14 12.118 11.286 13.047      NA
## West Baffin Bay                    NA     NA    NA     NA 13.857 11.136      NA
##                                   2018
## Amundsen Gulf Mouth                 NA
## Coronation Maud                15.1560
## Lancaster Sound                     NA
## Larsen Sound - Victoria Strait 10.7025
## Mackenzie Shelf                     NA
## NEG                                 NA
## North Water                    21.8750
## Peel Sound                          NA
## West Baffin Bay                     NA
cor(dfish$est_standard_length, dfish$open_water_day)
## [1] 0.4979243
plot(dfish$est_standard_length~dfish$open_water_day)

lattice::xyplot(est_standard_length ~ open_water_day|region, data = dfish)

#taille et carbone
plot(dfish$carbon_mg ~ dfish$fish_WW)

plot(dfish$carbon_mg ~ dfish$est_standard_length)

#taille et succès alimentaire
plot((feeding_success)~est_standard_length, data = dfish)

plot(log10(feeding_success)~est_standard_length, data = dfish)

summary(lm(log10(feeding_success) ~ est_standard_length, data = dfish))
## 
## Call:
## lm(formula = log10(feeding_success) ~ est_standard_length, data = dfish)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.41983 -0.26364  0.06344  0.31037  1.15673 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         -3.112901   0.058572 -53.147  < 2e-16 ***
## est_standard_length  0.009941   0.003028   3.283  0.00113 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4477 on 337 degrees of freedom
## Multiple R-squared:  0.031,  Adjusted R-squared:  0.02812 
## F-statistic: 10.78 on 1 and 337 DF,  p-value: 0.001134
plot(log10(feeding_success) ~ est_standard_length, data= dfish[dfish$est_standard_length>=25,]) #légère augmentation pour juvéniles en fonction de la condition, mais peu importe la condition chez les larves, le succès alimentaire reste assez variable (les résidus sont aussi moins dispersés)

plot(log10(feeding_success) ~ est_standard_length, data= dfish[dfish$est_standard_length<25,])

plot(log10(feeding_success) ~ fish_cond, data = dfish) #augmentation de FS

plot(log10(feeding_success) ~ fish_cond, data= dfish[dfish$est_standard_length>=25,])

plot(log10(feeding_success) ~ fish_cond, data= dfish[dfish$est_standard_length<25,]) 

plot(fish_cond~feeding_success, data = dfish)

plot(fish_cond~log10(feeding_success), data = dfish) # légère augmentation quand FS augmente

#pas vraiment de relation -> résidus pas assez dispersés (on n'a pas vraiment de valeurs extrêmes)

dfish$size_catego <- with(dfish, ifelse(est_standard_length<10, "<10", 
                                        ifelse(est_standard_length >=10 & est_standard_length<15, "10-15",
                                               ifelse(est_standard_length>=15 & est_standard_length<20, "15-20",
                                                      ifelse(est_standard_length>=20 & est_standard_length<25, "20-25",
                                                             ">25")))))
dfish$logFS <- log10(dfish$feeding_success)
xyplot(logFS ~ fish_cond|size_catego, data = dfish)

plot(fish_cond ~ est_standard_length, data = dfish) # la condition physique n'augmente pas vraiment avec la longueur du poisson

#FS et environnement
plot(logFS ~ NASC_zoo, data = dfish) #très légère augmentation mais rien de signifcatif

plot(logFS ~ surf_temp_degC, data = dfish) #légère augmentation

plot(logFS ~ surf_sal_kgm3, data = dfish) #pas grand chose

plot(logFS ~ open_water_day, data = dfish) #augmentation loglinéaire

dfish %>% select(logFS, open_water_day, surf_sal_kgm3, surf_temp_degC, NASC_zoo, fish_cond) %>% plot

#corrélation visible entre débâcle & température, débâcle et NASC_zoo et légère corrélation entre NASC_zoo et température
boxplot(logFS ~ region,data = dfish)

boxplot(logFS ~ province, data = dfish)

plot(feeding_success ~ region, data = dfish, las = 2) #pas mal similaire, mais toutes des valeurs bound à 0

plot(log10(feeding_success) ~ region, data = dfish, las = 2)

hist(dfish$feeding_success) #un beau neg binomial

hist(log10(dfish$feeding_success)) #normalish

lattice::bwplot(feeding_success ~ as.factor(year)|region, data = dfish)

lattice::bwplot(log10(feeding_success) ~ as.factor(year)|region, data = dfish) # variation LS, LV, MS, NW, PS

lattice::xyplot(log10(feeding_success) ~ (open_water_day)|region, data = dfish) #la variation individuelle est énorme toute région confondue

plot(feeding_success ~ as.factor(province), data = dfish)

plot(log10(feeding_success) ~ as.factor(province), data = dfish)

# Question 3 part 1

plot(prey_range_um/1000 ~ est_standard_length, data = dfish) # les grosses larves mangent à la fois de grosses et petites choses

plot(prey_min_um/1000 ~ est_standard_length, data = dfish)

plot(prey_max_um/1000 ~ est_standard_length, data = dfish)

plot(prey_median_um/1000 ~ est_standard_length, data = dfish) #vraiment rare que les larves, peu importe la taille, mangent de grosses affaires

plot(log10(feeding_success) ~ prey_median_um, data = dfish)

median(dfish[which(dfish$feeding_success > 0.001),"prey_median_um"])
## [1] 292.5
plot(log10(feeding_success) ~ prey_max_um, data = dfish)

plot(log10(feeding_success) ~ prey_min_um, data = dfish)

plot(log10(feeding_success) ~ prey_range_um, data = dfish)

plot((feeding_success) ~ prey_max_um, data = dfish)

plot((feeding_success) ~ prey_min_um, data = dfish)

plot((feeding_success) ~ prey_range_um, data = dfish)

#pas vraiment de relation entre taille des proies et FS

dfish$total_preys <- rowSums(dfish[,c(20:84)])
plot(total_preys ~ prey_max_um, data = dfish)

plot(total_preys ~ prey_min_um, data = dfish)

plot(total_preys~prey_range_um, data = dfish)

# Question 3 part 2

freqtable <- as.data.frame(table(gutcontent2$prey_id_comb,gutcontent2$region))
freqtable$Var2 <- ordered(freqtable$Var2, c("Mackenzie Shelf", "Amundsen Gulf Mouth", "Coronation Maud", "Larsen Sound - Victoria Strait", "Peel Sound", "Lancaster Sound", "North Water", "West Baffin Bay", "NEG"))
color = grDevices::colors()[grep('gr(a|e)y', grDevices::colors(), invert = T)]
colours_bar <- sample(color, size =15)
ggplot() +
  geom_bar(data =freqtable, mapping = aes(fill=Var1, y=Freq, x=Var2), position="fill", stat="identity") +
  scale_fill_manual(values =colours_bar)

freqtable2 <- as.data.frame(table(gutcontent2$prey_id_comb,gutcontent2$station))
ggplot() +
  geom_bar(data =freqtable2, mapping = aes(fill=Var1, y=Freq, x=Var2), position="fill", stat="identity") +
  scale_fill_manual(values =colours_bar)

carbon_input <- aggregate(carbon_mg ~ prey_id_comb + region, data = gutcontent2, FUN = sum)
ggplot() +
  geom_bar(data =carbon_input, mapping = aes(fill=prey_id_comb, y=carbon_mg, x=region), position="fill", stat="identity") +
  scale_fill_manual(values =colours_bar) + 
  labs(x = "", y = "")

carbon_input2 <- aggregate(carbon_mg ~ prey_id_comb + station, data = gutcontent2, FUN = sum)
ggplot() +
  geom_bar(data =carbon_input2, mapping = aes(fill=prey_id_comb, y=carbon_mg, x=station), position="fill", stat="identity") +
  scale_fill_manual(values =colours_bar) 

#en gros, bouffer des calanus = payant à souhait

PCA

fish_pca <-  dpreys_comb %>% select(-unique_fish_id,-fish_WW, -est_height_anus,-carbon_mg, -surf_sal_kgm3, -station, -sampling_date, -year, -bu_date, -c( 17:18, 20:42)) %>%PCA(quali.sup = c(2), quanti.sup = c(7))

dpreys_comb$pc1 <- fish_pca$ind$coord[, 1]
dpreys_comb$pc2 <- fish_pca$ind$coord[, 2]

pca.vars <- fish_pca$var$coord %>% data.frame
pca.vars$vars <- rownames(pca.vars)
pca.vars.m <- reshape2::melt(pca.vars, id.vars = "vars")
pca.var.sup <- fish_pca$quanti.sup$coord %>% data.frame

colours <- c("#6BCA54", "#FF5733", "#73C6B6", "#5E82D6", "#424949", "#A93226", "#C39BD3", "#FFC300", "#BDC3C7")
posArrowsX <- c(4, 4, 4, 3.25,4,4)
posArrowsY<- c(4.25, 3.5, 4.75, 4.25,4,4)
p_transf <- ggplot() +  
  geom_jitter(data = dpreys_comb, aes(x = pc1, y = pc2, colour = region), size = 2.5)+
  coord_equal()+
  theme_classic()+
  theme(panel.background = element_rect(fill = "transparent"), # bg of the panel
        plot.background = element_rect(fill = "transparent", color = NA),
        legend.position = c(0.25,0.19), 
        axis.text.x = element_text(size=20),
        axis.text.y = element_text(size = 20), 
        legend.text = element_text(size = 22), 
        axis.title = element_text(size = 25), 
        legend.title = element_text(size =22),
        legend.background = element_rect(fill=alpha('white', 0.4)))+
  labs(x= paste("PC1 (", round(fish_pca$eig["comp 1", "percentage of variance"], digits = 2), "%)"),
       y= paste("PC2 (", round(fish_pca$eig["comp 2", "percentage of variance"], digits = 2), "%)")) +
  geom_hline(yintercept = 0, lty = 2) +
  geom_vline(xintercept = 0, lty = 2) +
  geom_point(size = 2) +
  geom_segment(data = pca.vars, aes(x = 0, xend = Dim.1*4, y = 0, yend = Dim.2*4),
               arrow = arrow(length = unit(0.025, "npc"), type = "open"), 
               lwd = 0.5) + 
  geom_text(data = pca.vars, 
            aes(x = Dim.1*posArrowsX, y =  Dim.2*posArrowsY,
                label = c("LS","Jours depuis\ndébâcle","Abondance du\nzooplankton", "Température de\nsurface","Compé", "Condition physique")), 
            check_overlap = FALSE, size = 8.5)+ 
  scale_color_manual(values = colours, name = "Region") +
  geom_segment(data = pca.var.sup, aes(x = 0, xend = Dim.1*5, y = 0, yend = Dim.2*5),
             arrow = arrow(length = unit(0.025, "npc"), type = "open"),
             lwd = 0.5, lty = "dotted") + 
  geom_text(data = pca.var.sup, 
            aes(x = Dim.1*4, y =  Dim.2*4,
                label = "Succès\nalimentaire"), size = 8.5, col = '#868B8D')
p_transf
ggsave(file="/Users/pascalecaissy/Dropbox/MSc 2019-2020/bosaFS/PCA_fish.png", plot = p_transf, dpi = 500, width = 14, height = 10, bg = "transparent")

CA

CA.preys <- dpreys_comb %>% select(22:42) %>% CA()
summary(CA.preys) #first two axes explains 24%, not really good
plot(CA.preys)

# PCA.prey <- dpreys_comb %>% select(-unique_fish_id, -est_height_anus, -fish_WW, -carbon_mg) %>% PCA(quanti.sup = c(7), quali.sup = c(2))

env <- dpreys_comb %>% select(-region, -unique_fish_id, -est_height_anus, -fish_WW, -carbon_mg, -surf_sal_kgm3, -est_standard_length, -c(13:34))
preys <- dpreys_comb %>% select(13:34) %>% decostand(method = "hellinger")
rda.preys <- rda(X = preys, Y = env, scale = TRUE)
plot(rda.preys)
summary(rda.preys) #not really good, two first axes explain only 11%...

preys_pres_abs <-dpreys_comb %>% select(13:34) %>% sapply(FUN = function(x) 
  ifelse(x > 0, 1, 0))
rda.preys.presabs <- rda(X = preys_pres_abs, Y = env, scale = TRUE)
plot(rda.preys.presabs)
summary(rda.preys.presabs)

PCA.prey <- dpreys_comb %>% select(-unique_fish_id, -est_height_anus, -fish_WW, -carbon_mg, -est_standard_length) %>% PCA(quanti.sup = c(6), quali.sup = c(1))

Bray Curtis, MarineCusa, and others…….

wo data sets were generated for the multivariate analysis: Dataset A, where the gravimetric contribution of BDC in percent (%W) was calculated for each stomach by dividing prey taxa weight by total prey weight for each stomach and multiplying the result by 100;

and dataset B, in which the average BDC weight was calculated for each trawl by adding total prey taxa weight for each trawl, dividing total prey taxa weight by total prey weight for each trawl, and multiplying the result by 100. Dataset B calculation were performed separately for each polar cod size category

Sélection de modèles

#vérification ajustement des modèles lm----
lm.fs.untranf <- lm(feeding_success ~ open_water_day + fish_cond, data = dfish, na.action = na.omit)
plot(lm.fs.untranf)

lm.glob1 <- lm(log10(feeding_success) ~ open_water_day + fish_cond, data = dfish, na.action = na.omit)
summary(lm.glob1)
## 
## Call:
## lm(formula = log10(feeding_success) ~ open_water_day + fish_cond, 
##     data = dfish, na.action = na.omit)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.45771 -0.25781  0.05526  0.28660  1.24458 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -3.115499   0.041676 -74.755  < 2e-16 ***
## open_water_day  0.003567   0.000692   5.155 4.33e-07 ***
## fish_cond       0.183561   0.073807   2.487   0.0134 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4321 on 336 degrees of freedom
## Multiple R-squared:  0.1003, Adjusted R-squared:  0.09497 
## F-statistic: 18.74 on 2 and 336 DF,  p-value: 1.932e-08
#pdf("/Users/pascalecaissy/Dropbox/MSc 2019-2020/bosaFS/glob1_res.pdf")
plot(rstudent(lm.glob1) ~ fitted(lm.glob1), ylab = "Résidus de Student", xlab = "Valeurs prédites")

#dev.off()

plot(hatvalues(lm.glob1), 
     ylab = "Hat values",
     xlab = "Observations", main = "Effet de levier")
abline(h = 2*3/324, lty = 2) #on a un effet de levier pour certaines valeurs, mais je ne crois pas qu'il soit trop

# pdf("/Users/pascalecaissy/Dropbox/MSc 2019-2020/bosaFS/glob1_cook.pdf")
plot(cooks.distance(lm.glob1),
     ylab = "Distance de Cook", 
     xlab = "Observations", 
     main = "Influence des observations")
abline(h = 3/(324-4), lty = 2) #rien en haut de 1, mais un bon nombre sont en haut du seuil

# dev.off()

lm.glob2 <- lm(log10(feeding_success) ~ surf_temp_degC + NASC_zoo + fish_cond, data = dfish, na.action = na.omit)
summary(lm.glob2)
## 
## Call:
## lm(formula = log10(feeding_success) ~ surf_temp_degC + NASC_zoo + 
##     fish_cond, data = dfish, na.action = na.omit)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.4832 -0.2482  0.0573  0.3056  1.1629 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -2.9486490  0.0407859 -72.296  < 2e-16 ***
## surf_temp_degC  0.0024438  0.0132247   0.185  0.85350    
## NASC_zoo        0.0004823  0.0039241   0.123  0.90226    
## fish_cond       0.2391899  0.0791796   3.021  0.00271 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4494 on 335 degrees of freedom
## Multiple R-squared:  0.02951,    Adjusted R-squared:  0.02082 
## F-statistic: 3.395 on 3 and 335 DF,  p-value: 0.01815
# pdf("/Users/pascalecaissy/Dropbox/MSc 2019-2020/bosaFS/glob2_res.pdf")
plot(rstudent(lm.glob2) ~ fitted(lm.glob2), ylab = "Résidus de Student", xlab = "Valeurs prédites")

# dev.off()

plot(hatvalues(lm.glob2), ylab = "Hat values",
     xlab = "Observations", main = "Effet de levier")
abline(h = 2*4/324, lty = 2) #on a un effet de levier pour certaines valeurs, mais je ne crois pas qu'il soit trop

# pdf("/Users/pascalecaissy/Dropbox/MSc 2019-2020/bosaFS/glob2_cook.pdf")
plot(cooks.distance(lm.glob2), ylab = "Distance de Cook", xlab = "Observations", main = "Influence des observations")
abline(h = 4/(324-4), lty = 2) #rien en haut de 1, mais un bon nombre sont en haut du seuil

dev.off()
## null device 
##           1
#modèles mixtes----
lme.null <- lme(log10(feeding_success) ~1, random = ~1|region, data = dfish, na.action = na.omit, method = "ML")
summary(lme.null)
## Linear mixed-effects model fit by maximum likelihood
##  Data: dfish 
##        AIC      BIC    logLik
##   412.5423 424.0203 -203.2712
## 
## Random effects:
##  Formula: ~1 | region
##         (Intercept)  Residual
## StdDev:   0.1415781 0.4318548
## 
## Fixed effects: log10(feeding_success) ~ 1 
##                 Value  Std.Error  DF   t-value p-value
## (Intercept) -2.913366 0.05418298 330 -53.76903       0
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.2166700 -0.5841590  0.1201451  0.6732820  3.0226108 
## 
## Number of Observations: 339
## Number of Groups: 9
lme.glob1 <- lme(log10(feeding_success) ~ surf_temp_degC + NASC_zoo + fish_cond + indstandard, random = ~1|region, data = dfish, na.action = na.omit, method = "ML")
summary(lme.glob1) #serait mieux d'utiliser une transformation log10
## Linear mixed-effects model fit by maximum likelihood
##  Data: dfish 
##        AIC      BIC    logLik
##   413.8603 440.6423 -199.9302
## 
## Random effects:
##  Formula: ~1 | region
##         (Intercept)  Residual
## StdDev:   0.1347363 0.4279663
## 
## Fixed effects: log10(feeding_success) ~ surf_temp_degC + NASC_zoo + fish_cond +      indstandard 
##                     Value  Std.Error  DF   t-value p-value
## (Intercept)    -2.8890499 0.06751992 326 -42.78811  0.0000
## surf_temp_degC  0.0008893 0.01568588 326   0.05670  0.9548
## NASC_zoo       -0.0036067 0.00477970 326  -0.75458  0.4510
## fish_cond       0.1759375 0.07905116 326   2.22562  0.0267
## indstandard    -0.1204569 0.22282177 326  -0.54060  0.5892
##  Correlation: 
##                (Intr) srf__C NASC_z fsh_cn
## surf_temp_degC -0.480                     
## NASC_zoo       -0.017 -0.588              
## fish_cond       0.044 -0.179  0.191       
## indstandard    -0.276  0.132  0.007 -0.021
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.19707455 -0.51683825  0.09907741  0.64911501  3.07875340 
## 
## Number of Observations: 339
## Number of Groups: 9
lme.glob2 <- lme(log10(feeding_success) ~ open_water_day + fish_cond, random = ~1|region, data = dfish, na.action = na.omit, method = "ML")
summary(lme.glob2)
## Linear mixed-effects model fit by maximum likelihood
##  Data: dfish 
##        AIC      BIC    logLik
##   392.4367 411.5667 -191.2183
## 
## Random effects:
##  Formula: ~1 | region
##         (Intercept)  Residual
## StdDev:   0.1226326 0.4176763
## 
## Fixed effects: log10(feeding_success) ~ open_water_day + fish_cond 
##                    Value  Std.Error  DF   t-value p-value
## (Intercept)    -3.165830 0.07516693 328 -42.11732  0.0000
## open_water_day  0.004910 0.00113264 328   4.33509  0.0000
## fish_cond       0.153737 0.07540482 328   2.03882  0.0423
##  Correlation: 
##                (Intr) opn_w_
## open_water_day -0.766       
## fish_cond       0.070 -0.097
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.1276007 -0.5938883  0.1068560  0.6323493  3.0566429 
## 
## Number of Observations: 339
## Number of Groups: 9
performance::r2(lme.glob2)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.217
##      Marginal R2: 0.150
lme.glob3 <- lme(log10(feeding_success) ~ open_water_day + fish_cond +surf_temp_degC + NASC_zoo + indstandard, random = ~1|region, data = dfish, na.action = na.omit, method = "ML")
summary(lme.glob3)
## Linear mixed-effects model fit by maximum likelihood
##  Data: dfish 
##        AIC      BIC    logLik
##   392.1609 422.7689 -188.0805
## 
## Random effects:
##  Formula: ~1 | region
##         (Intercept)  Residual
## StdDev:   0.1486085 0.4121206
## 
## Fixed effects: log10(feeding_success) ~ open_water_day + fish_cond + surf_temp_degC +      NASC_zoo + indstandard 
##                    Value  Std.Error  DF   t-value p-value
## (Intercept)    -3.182257 0.09219465 325 -34.51673  0.0000
## open_water_day  0.006638 0.00134656 325   4.92936  0.0000
## fish_cond       0.136769 0.07680118 325   1.78082  0.0759
## surf_temp_degC -0.013225 0.01558379 325  -0.84865  0.3967
## NASC_zoo       -0.006531 0.00471398 325  -1.38551  0.1668
## indstandard     0.174968 0.22380914 325   0.78177  0.4349
##  Correlation: 
##                (Intr) opn_w_ fsh_cn srf__C NASC_z
## open_water_day -0.653                            
## fish_cond       0.090 -0.096                     
## surf_temp_degC -0.217 -0.183 -0.149              
## NASC_zoo        0.063 -0.114  0.194 -0.553       
## indstandard    -0.358  0.264 -0.045  0.067 -0.015
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.21435875 -0.60432496  0.08542066  0.64163558  3.16442861 
## 
## Number of Observations: 339
## Number of Groups: 9
performance::r2(lme.glob3)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.286
##      Marginal R2: 0.193
lme.bu <- lme(log10(feeding_success) ~ open_water_day, random = ~1|region, data = dfish, na.action = na.omit, method = "ML")
summary(lme.bu)
## Linear mixed-effects model fit by maximum likelihood
##  Data: dfish 
##        AIC      BIC    logLik
##   394.5949 409.8989 -193.2975
## 
## Random effects:
##  Formula: ~1 | region
##         (Intercept)  Residual
## StdDev:   0.1280725 0.4199369
## 
## Fixed effects: log10(feeding_success) ~ open_water_day 
##                    Value  Std.Error  DF   t-value p-value
## (Intercept)    -3.179778 0.07694985 329 -41.32274       0
## open_water_day  0.005202 0.00114956 329   4.52555       0
##  Correlation: 
##                (Intr)
## open_water_day -0.761
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.16166520 -0.57445494  0.09310238  0.66866541  2.98507700 
## 
## Number of Observations: 339
## Number of Groups: 9
lme.cond <- lme(log10(feeding_success) ~ fish_cond, random = ~1|region, data = dfish, na.action = na.omit, method = "ML")
summary(lme.cond)
## Linear mixed-effects model fit by maximum likelihood
##  Data: dfish 
##        AIC      BIC    logLik
##   408.8741 424.1781 -200.4371
## 
## Random effects:
##  Formula: ~1 | region
##         (Intercept)  Residual
## StdDev:   0.1305161 0.4288933
## 
## Fixed effects: log10(feeding_success) ~ fish_cond 
##                  Value  Std.Error  DF   t-value p-value
## (Intercept) -2.9155946 0.05090030 329 -57.28050  0.0000
## fish_cond    0.1845831 0.07704973 329   2.39564  0.0171
##  Correlation: 
##           (Intr)
## fish_cond -0.006
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.16832072 -0.54686777  0.09007215  0.65614797  3.08828777 
## 
## Number of Observations: 339
## Number of Groups: 9
lme.comp <- lme(log10(feeding_success) ~ indstandard, random = ~1|region, data = dfish, na.action = na.omit, method = "ML")
summary(lme.comp)
## Linear mixed-effects model fit by maximum likelihood
##  Data: dfish 
##        AIC      BIC    logLik
##   414.3573 429.6613 -203.1786
## 
## Random effects:
##  Formula: ~1 | region
##         (Intercept)  Residual
## StdDev:   0.1399442 0.4318376
## 
## Fixed effects: log10(feeding_success) ~ indstandard 
##                  Value  Std.Error  DF   t-value p-value
## (Intercept) -2.9084324 0.05509209 329 -52.79220  0.0000
## indstandard -0.0951137 0.22091146 329  -0.43055  0.6671
##  Correlation: 
##             (Intr)
## indstandard -0.217
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.2231890 -0.5886669  0.1248675  0.6830543  2.9948472 
## 
## Number of Observations: 339
## Number of Groups: 9
lme.NASC <- lme(log10(feeding_success) ~ NASC_zoo, random = ~1|region, data = dfish, na.action = na.omit, method = "ML")
summary(lme.NASC)
## Linear mixed-effects model fit by maximum likelihood
##  Data: dfish 
##        AIC      BIC    logLik
##   413.3037 428.6077 -202.6519
## 
## Random effects:
##  Formula: ~1 | region
##         (Intercept)  Residual
## StdDev:   0.1474387 0.4306864
## 
## Fixed effects: log10(feeding_success) ~ NASC_zoo 
##                  Value  Std.Error  DF   t-value p-value
## (Intercept) -2.8862064 0.06073509 329 -47.52124  0.0000
## NASC_zoo    -0.0043282 0.00386754 329  -1.11910  0.2639
##  Correlation: 
##          (Intr)
## NASC_zoo -0.388
## 
## Standardized Within-Group Residuals:
##         Min          Q1         Med          Q3         Max 
## -3.24785165 -0.57784643  0.08928769  0.65987904  3.05791490 
## 
## Number of Observations: 339
## Number of Groups: 9
lme.temp <- lme(log10(feeding_success) ~ surf_temp_degC, random = ~1|region, data = dfish, na.action = na.omit, method = "ML")
summary(lme.temp)
## Linear mixed-effects model fit by maximum likelihood
##  Data: dfish 
##        AIC      BIC   logLik
##   414.4939 429.7979 -203.247
## 
## Random effects:
##  Formula: ~1 | region
##         (Intercept)  Residual
## StdDev:   0.1430711 0.4317298
## 
## Fixed effects: log10(feeding_success) ~ surf_temp_degC 
##                     Value  Std.Error  DF   t-value p-value
## (Intercept)    -2.9045914 0.06695295 329 -43.38258  0.0000
## surf_temp_degC -0.0027957 0.01259061 329  -0.22205  0.8244
##  Correlation: 
##                (Intr)
## surf_temp_degC -0.577
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.2347463 -0.5854568  0.1171700  0.6707042  3.0315922 
## 
## Number of Observations: 339
## Number of Groups: 9
list.models <- list(null = lme.null, tempNASCFishcond = lme.glob1, buFishcond= lme.glob2, buOnly = lme.bu, tempONLY = lme.temp, condOnly = lme.cond, NASConly = lme.NASC, compONLY = lme.comp)

selection <- AICcmodavg::aictab(list.models)
## Warning: no function found corresponding to methods exports from 'raster' for:
## 'wkt'
selection #les deux premiers modèles semblent être ceux qui ont le plus de poids
## 
## Model selection based on AICc:
## 
##                  K   AICc Delta_AICc AICcWt Cum.Wt      LL
## buFishcond       5 392.62       0.00   0.74   0.74 -191.22
## buOnly           4 394.71       2.10   0.26   1.00 -193.30
## condOnly         4 408.99      16.38   0.00   1.00 -200.44
## null             3 412.61      20.00   0.00   1.00 -203.27
## NASConly         4 413.42      20.81   0.00   1.00 -202.65
## tempNASCFishcond 7 414.20      21.58   0.00   1.00 -199.93
## compONLY         4 414.48      21.86   0.00   1.00 -203.18
## tempONLY         4 414.61      22.00   0.00   1.00 -203.25
est.bu <- AICcmodavg::modavg(list.models, "open_water_day")
est.bu
## 
## Multimodel inference on "open_water_day" based on AICc
## 
## AICc table used to obtain model-averaged estimate:
## 
##            K   AICc Delta_AICc AICcWt Estimate SE
## buFishcond 5 392.62        0.0   0.74     0.00  0
## buOnly     4 394.71        2.1   0.26     0.01  0
## 
## Model-averaged estimate: 0 
## Unconditional SE: 0 
## 95% Unconditional confidence interval: 0, 0.01
est.fish.cond <- AICcmodavg::modavg(list.models, "fish_cond")
est.fish.cond
## 
## Multimodel inference on "fish_cond" based on AICc
## 
## AICc table used to obtain model-averaged estimate:
## 
##                  K   AICc Delta_AICc AICcWt Estimate   SE
## tempNASCFishcond 7 414.20      21.58      0     0.18 0.08
## buFishcond       5 392.62       0.00      1     0.15 0.08
## condOnly         4 408.99      16.38      0     0.18 0.08
## 
## Model-averaged estimate: 0.15 
## Unconditional SE: 0.08 
## 95% Unconditional confidence interval: 0.01, 0.3
est.temp <- AICcmodavg::modavgShrink(list.models, "surf_temp_degC")
## Warning in modavgShrink.AIClme(list.models, "surf_temp_degC"): 
## Variables do not appear with same frequency across models, proceed with caution
est.temp
## 
## Multimodel inference on "surf_temp_degC" based on AICc
## 
## AICc table used to obtain model-averaged estimate with shrinkage:
## 
##                  K   AICc Delta_AICc AICcWt Estimate   SE
## null             3 412.61      20.00   0.00        0 0.00
## tempNASCFishcond 7 414.20      21.58   0.00        0 0.02
## buFishcond       5 392.62       0.00   0.74        0 0.00
## buOnly           4 394.71       2.10   0.26        0 0.00
## tempONLY         4 414.61      22.00   0.00        0 0.01
## condOnly         4 408.99      16.38   0.00        0 0.00
## NASConly         4 413.42      20.81   0.00        0 0.00
## compONLY         4 414.48      21.86   0.00        0 0.00
## 
## Model-averaged estimate with shrinkage: 0 
## Unconditional SE: 0 
## 95% Unconditional confidence interval: 0, 0
est.zoo <- AICcmodavg::modavgShrink(list.models, "NASC_zoo")
## Warning in modavgShrink.AIClme(list.models, "NASC_zoo"): 
## Variables do not appear with same frequency across models, proceed with caution
est.zoo
## 
## Multimodel inference on "NASC_zoo" based on AICc
## 
## AICc table used to obtain model-averaged estimate with shrinkage:
## 
##                  K   AICc Delta_AICc AICcWt Estimate SE
## null             3 412.61      20.00   0.00        0  0
## tempNASCFishcond 7 414.20      21.58   0.00        0  0
## buFishcond       5 392.62       0.00   0.74        0  0
## buOnly           4 394.71       2.10   0.26        0  0
## tempONLY         4 414.61      22.00   0.00        0  0
## condOnly         4 408.99      16.38   0.00        0  0
## NASConly         4 413.42      20.81   0.00        0  0
## compONLY         4 414.48      21.86   0.00        0  0
## 
## Model-averaged estimate with shrinkage: 0 
## Unconditional SE: 0 
## 95% Unconditional confidence interval: 0, 0
new.data <- data.frame(
  open_water_day = seq(from = min(dfish$open_water_day), to = max(dfish$open_water_day), length.out = 100),
  fish_cond = mean(dfish$fish_cond),
  region = "Amundsen Gulf Mouth",
  surf_temp_degC = mean(dfish$surf_temp_degC),
  NASC_zoo = mean(dfish$NASC_zoo),
  indstandard = mean(dfish$indstandard)
)
preds <- AICcmodavg::modavgPred(cand.set = list.models, newdata = new.data)
preds2 <- predict(lme.glob2, newdata =new.data)

new.data$fit <- preds$mod.avg.pred
new.data$se <- preds$uncond.se
new.data$low95 <- preds$lower.CL
new.data$supp95 <- preds$upper.CL

#pdf("/Users/pascalecaissy/Dropbox/MSc 2019-2020/bosaFS/preds_bu.pdf")
plot(smooth(fit) ~ open_water_day, data = new.data,
     ylab = "log10(Succès alimentaire des jeunes morues polaires)",
     xlab = "Nombre de jours juliens depuis la débâcle",
     main = "Valeurs prédites pondérées de succès alimentaire (+- IC 95%)",
     type = 'l')
lines(x = new.data$open_water_day, y = smooth(new.data$supp95), lty = 'dotted')
lines(x = new.data$open_water_day, y = smooth(new.data$low95), lty = 'dotted')
#dev.off()

Poster RSA QO

mean_logFS <- aggregate(logFS~ region + open_water_day, data = dfish, FUN = mean)
mean_logFS$region <- ordered(mean_logFS$region, c("Mackenzie Shelf", "Amundsen Gulf Mouth", "Coronation Maud", "Larsen Sound - Victoria Strait", "Peel Sound", "Lancaster Sound", "North Water", "West Baffin Bay", "NEG"))

colours <- pals::stepped(n = 24)[c(1,4,9,11,12,13, 15,16,17)]
FS_OWD <- ggplot() +
  geom_point(data = mean_logFS, mapping = aes(x = open_water_day, y = logFS, color = region), size = 10)+
  geom_line(data = new.data, mapping = aes(y = fit, x = open_water_day)) + 
  geom_ribbon(data = new.data, aes(x = open_water_day, y = NULL, ymin = low95, ymax = supp95),
                 alpha = .15) + 
  theme_classic()+
  theme(panel.background = element_rect(fill = "transparent"), # bg of the panel
        plot.background = element_rect(fill = "transparent", color = NA),
        legend.position = c(0.25,0.78), 
        axis.text.x = element_text(size=30),
        axis.text.y = element_text(size = 30), 
        legend.text = element_text(size = 22), 
        axis.title = element_text(size = 30), 
        legend.title = element_text(size =22),
        legend.background = element_rect(fill="transparent")) + 
  labs(x = "Nombre de jours depuis la débâcle", y = "log10(Succès alimentaire moyen)") + 
  scale_color_manual(values = colours, name = "Région", labels = c("Mackenzie Shelf", "Amundsen Gulf Mouth", "Coronation Maud", "Larsen Sound - Victoria Strait", "Peel Sound", "Lancaster Sound", "North Water", "West Baffin Bay", "North-Eastern\nGreenland Sea")) + 
  scale_x_continuous(breaks = c(-20, 0, 20, 40, 60, 80, 100))

FS_OWD
ggsave(file="/Users/pascalecaissy/Dropbox/MSc 2019-2020/bosaFS/FS_OWD_fish.png", plot = FS_OWD, dpi = 500, width = 14, height = 10, bg = "transparent")
library(pals)
colours_bar <- pals::stepped(n = 24)[c(20:17, 16:13, 12,5,8,4:1)]

gutcontent2$OWD_cat <- with(gutcontent2, ifelse(open_water_day <= 0, "before_bu",
                                              ifelse(open_water_day >0 & open_water_day <20, "0-20",
                                                     ifelse(open_water_day >=20 & open_water_day<40, "20-40",
                                                            ifelse(open_water_day >=40 & open_water_day <60, "40-60",
                                                                   ifelse(open_water_day >=60 & open_water_day < 80, "60-80", "80+"))))))

freqtable <- as.data.frame(table(gutcontent2$prey_id_comb,gutcontent2$region))
freqtable$Var2 <- ordered(freqtable$Var2, c("Mackenzie Shelf", "Amundsen Gulf Mouth", "Coronation Maud", "Larsen Sound - Victoria Strait", "Peel Sound", "Lancaster Sound", "North Water", "West Baffin Bay", "NEG"))
freqtable$Var1 <- ordered(freqtable$Var1, c("appendicularia", "bivalve larvae", "egg", "limacina helicina", "calanoid sp c", "cyclopoidae sp c",  "other copepodite sp", "other copepodite n", "other", "pseudocalanus sp c", "pseudocalanus sp n", "calanus glacialis c", "calanus hyperboreus c", "calanus sp c", "calanus sp n"))


preys_occ <- ggplot() +
  geom_bar(data =freqtable, mapping = aes(fill=Var1, y=Freq, x=Var2), position="fill", stat="identity") +
  ggtitle( "Fréquence d'occurence\ndans la diète")+
   scale_fill_manual(name = "", values =colours_bar, labels = c("Appendicularia", 
                                                    "Bivalvia", 
                                                    "Oeuf", 
                                                    "Gastropoda", 
                                                    "Calanoidae spp.", 
                                                    "Cyclopoidae spp.", 
                                                    "Autre copépodite", 
                                                    "Autre nauplii",
                                                    "Autre type de proies",
                                                    expression(paste(italic("Pseudocalanus "), "spp.", sep = "")),
                                                    expression(paste(italic("Pseudocalanus "), "spp.  Nauplii", sep = "")),
                                                    expression(italic("Calanus glacialis")), 
                                                    expression(italic("Calanus hyperboreus")), 
                                                    expression(paste(italic("Calanus "), "spp.", sep = "")), 
                                                    expression(paste(italic("Calanus "), "spp.  Nauplii", sep = ""))
                                                    ))+ 
  theme_classic() +
  theme(panel.background = element_rect(fill = "transparent"), # bg of the panel
        plot.background = element_rect(fill = "transparent", color = NA),
        axis.text.x = element_text(size=20),
        axis.text.y = element_text(size = 20), 
        legend.text = element_text(size = 22),
        legend.text.align = 0,
        axis.title = element_text(size = 25), 
        legend.title = element_text(size =22),
        legend.background = element_rect(fill=alpha('white', 0.4)),
        plot.title = element_text(size = 30, hjust = 0.5)) + 
  labs(y = "Proportion", x = "") + 
  scale_x_discrete(labels = c("MS", "AGM", "CM", "LV", "PS", "LS", "NW", "WBB", "NEG"))
preys_occ
#ggsave(file="/Users/pascalecaissy/Dropbox/MSc 2019-2020/bosaFS/preys_occ.png", plot = preys_occ, dpi = 500, width = 14, height = 10, bg = "transparent")

carbon_input <- aggregate(carbon_mg ~ prey_id_comb + region, data = gutcontent2, FUN = sum)
carbon_input$region <- ordered(carbon_input$region, c("Mackenzie Shelf", "Amundsen Gulf Mouth", "Coronation Maud", "Larsen Sound - Victoria Strait", "Peel Sound", "Lancaster Sound", "North Water", "West Baffin Bay", "NEG"))
carbon_input$prey_id_comb <- ordered(carbon_input$prey_id_comb, c("appendicularia", "bivalve larvae", "egg", "limacina helicina", "calanoid sp c", "cyclopoidae sp c",  "other copepodite sp", "other copepodite n", "other", "pseudocalanus sp c", "pseudocalanus sp n", "calanus glacialis c", "calanus hyperboreus c", "calanus sp c", "calanus sp n"))
preys_carbon <- ggplot() +
  geom_bar(data =carbon_input, mapping = aes(fill=prey_id_comb, y=carbon_mg, x=region), position="fill", stat="identity") +
  ggtitle("Proportion en carbone\ndans la diète")+
   scale_fill_manual(name = "", values =colours_bar, labels = c("Appendicularia", 
                                                    "Bivalvia", 
                                                    "Oeuf", 
                                                    "Gastropoda", 
                                                    "Calanoidae spp.", 
                                                    "Cyclopoidae spp.", 
                                                    "Autre copépodite", 
                                                    "Autre nauplii",
                                                    "Autre type de proies",
                                                    expression(paste(italic("Pseudocalanus "), "spp.", sep = "")),
                                                    expression(paste(italic("Pseudocalanus "), "spp.  Nauplii", sep = "")),
                                                    expression(italic("Calanus glacialis")), 
                                                    expression(italic("Calanus hyperboreus")), 
                                                    expression(paste(italic("Calanus "), "spp.", sep = "")), 
                                                    expression(paste(italic("Calanus "), "spp.  Nauplii", sep = ""))
                                                    )) + 
  theme_classic() +
  theme(panel.background = element_rect(fill = "transparent"), # bg of the panel
        plot.background = element_rect(fill = "transparent", color = NA),
        axis.text.x = element_text(size=20),
        axis.text.y = element_text(size = 20), 
        legend.text = element_text(size = 22),
        legend.text.align = 0,
        axis.title = element_text(size = 25), 
        legend.title = element_text(size =22),
        legend.background = element_rect(fill=alpha('white', 0.4)),
        plot.title = element_text(size = 30, hjust = 0.5)) + 
  labs(y = "", x = "Région") + 
  scale_x_discrete(labels = c("MS", "AGM", "CM", "LV", "PS", "LS", "NW", "WBB", "NEG"))
preys_carbon

#ggsave(file="/Users/pascalecaissy/Dropbox/MSc 2019-2020/bosaFS/preys_carbon.png", plot = preys_carbon, dpi = 500, width = 14, height = 10, bg = "transparent")



freqtable3 <- as.data.frame(table(gutcontent2$prey_id_comb,gutcontent2$OWD_cat))
freqtable3$Var1 <- ordered(freqtable3$Var1, c("appendicularia", "bivalve larvae", "egg", "limacina helicina", "calanoid sp c", "cyclopoidae sp c",  "other copepodite sp", "other copepodite n", "other", "pseudocalanus sp c", "pseudocalanus sp n", "calanus glacialis c", "calanus hyperboreus c", "calanus sp c", "calanus sp n"))
freqtable3$Var2 <- ordered(freqtable3$Var2,c("before_bu", "0-20","20-40", "40-60", "60-80", "80+"))


preys_occ_owd <- ggplot() +
  geom_bar(data =freqtable3, mapping = aes(fill=Var1, y=Freq, x=Var2), position="fill", stat="identity") +
  ggtitle( "Fréquence d'occurence\ndes proies dans la diète")+
   scale_fill_manual(name = "", values =colours_bar, labels = c("Appendicularia", 
                                                    "Bivalvia", 
                                                    "Oeuf", 
                                                    "Gastropoda", 
                                                    "Calanoidae spp.", 
                                                    "Cyclopoidae spp.", 
                                                    "Autre copépodite", 
                                                    "Autre nauplii",
                                                    "Autre type de proies",
                                                    expression(paste(italic("Pseudocalanus "), "spp.", sep = "")),
                                                    expression(paste(italic("Pseudocalanus "), "spp. Nauplii", sep = "")),
                                                    expression(italic("Calanus glacialis")), 
                                                    expression(italic("Calanus hyperboreus")), 
                                                    expression(paste(italic("Calanus "), "spp.", sep = "")), 
                                                    expression(paste(italic("Calanus "), "spp. Nauplii", sep = ""))
                                                    ))+ 
  theme_classic() +
  theme(panel.background = element_rect(fill = "transparent"), # bg of the panel
        plot.background = element_rect(fill = "transparent", color = NA),
        axis.text.x = element_text(size=30),
        axis.text.y = element_text(size = 30), 
        legend.text = element_text(size = 22),
        legend.text.align = 0,
        axis.title = element_text(size = 35), 
        legend.title = element_text(size =22),
        legend.background = element_rect(fill=alpha('white', 0.4)),
        plot.title = element_text(size = 35, hjust = 0.5)) + 
  labs(y = "Proportion", x = "Nombre de jours depuis la débâcle")+
  scale_x_discrete(labels = c("Avant\ndébâcle", "0-20 j.", "20-40 j.", "40-60 j.", "60-80 j.", "80+ j."))
preys_occ_owd
ggsave(file="/Users/pascalecaissy/Dropbox/MSc 2019-2020/bosaFS/preys_OWD.png", plot = preys_occ_owd, dpi = 500, width = 17, height = 14, bg = "transparent")

carbon_input2 <- aggregate(carbon_mg ~ prey_id_comb + OWD_cat + unique_fish_id, data = gutcontent2, FUN = sum)
carbon_input3 <- aggregate(carbon_mg ~ prey_id_comb + OWD_cat, data = carbon_input2, FUN=mean)
carbon_input3$OWD_cat <- ordered(carbon_input3$OWD_cat ,c("before_bu", "0-20","20-40", "40-60", "60-80", "80+"))
carbon_input3$prey_id_comb <- ordered(carbon_input3$prey_id_comb, c("appendicularia", "bivalve larvae", "egg", "limacina helicina", "calanoid sp c", "cyclopoidae sp c",  "other copepodite sp", "other copepodite n", "other", "pseudocalanus sp c", "pseudocalanus sp n", "calanus glacialis c", "calanus hyperboreus c", "calanus sp c", "calanus sp n"))
preys_carbon_OWD <- ggplot() +
  geom_bar(data =carbon_input3, mapping = aes(y=carbon_mg, x=OWD_cat, fill = prey_id_comb), stat="identity") +
  ggtitle("Moyenne de carbone ingéré par poisson\net contribution des proies")+
  scale_fill_manual(name = "", values =colours_bar, labels = c("Appendicularia", 
                                                    "Bivalvia", 
                                                    "Oeuf", 
                                                    "Gastropoda", 
                                                    "Calanoidae spp.", 
                                                    "Cyclopoidae spp.", 
                                                    "Autre copépodite", 
                                                    "Autre nauplii",
                                                    "Autre type de proies",
                                                    expression(paste(italic("Pseudocalanus "), "spp.", sep = "")),
                                                    expression(paste(italic("Pseudocalanus "), "spp. Nauplii", sep = "")),
                                                    expression(italic("Calanus glacialis")), 
                                                    expression(italic("Calanus hyperboreus")), 
                                                    expression(paste(italic("Calanus "), "spp.", sep = "")), 
                                                    expression(paste(italic("Calanus "), "spp. Nauplii", sep = ""))
                                                    ))+ 
  theme_classic() +
  theme(panel.background = element_rect(fill = "transparent"), # bg of the panel
        plot.background = element_rect(fill = "transparent", color = NA),
        axis.text.x = element_text(size=28),
        axis.text.y = element_text(size = 28), 
        legend.text = element_text(size = 25),
        legend.text.align = 0,
        axis.title = element_text(size = 35), 
        legend.title = element_text(size =22),
        legend.background = element_rect(fill=alpha('white', 0.4)),
        plot.title = element_text(size = 35, hjust = 0.5))+
  scale_x_discrete(labels = c("Avant\ndébâcle", "0-20 j.", "20-40 j.", "40-60 j.", "60-80 j.", "80+"))+
  labs(y = "Moyenne de carbone ingéré/poisson (mg)", x = "Nombre de jours depuis la débâcle")
preys_carbon_OWD
ggsave(file="/Users/pascalecaissy/Dropbox/MSc 2019-2020/bosaFS/preys_carbon_OWD.png", plot = preys_carbon_OWD, dpi = 500, width = 17, height = 14, bg = "transparent")
#le faire pour des catégories de succès alimentaire


carbon_input4 <- aggregate(carbon_mg ~ prey_id_comb + open_water_day + unique_fish_id, data = gutcontent2, FUN = sum)
carbon_input5 <- aggregate(carbon_mg ~ prey_id_comb + open_water_day, data = carbon_input4, FUN = mean)
carbon_input6 <- aggregate(carbon_mg ~ open_water_day, data = carbon_input5, FUN = sum)
carbon_input_caspC <- aggregate(carbon_mg ~ open_water_day, 
                                data = carbon_input4[carbon_input5$prey_id_comb %in% 
                                                       c("calanus glacialis c", 
                                                         "calanus hyperboreus c", 
                                                         "calanus sp c"),], FUN = sum)

ggplot()+
  geom_point(data = carbon_input6, mapping = aes(x=open_water_day, y = carbon_mg), color = "black")+
  geom_point(data = carbon_input_caspC, mapping = aes(x=open_water_day, y = carbon_mg), color = "#990F26")+
  geom_smooth(method = "loess", span = 0.75,data = carbon_input6, mapping = aes(x=open_water_day, y = carbon_mg), color = "black", fill = "grey") +
  geom_smooth(method = "loess", span = 0.75, data = carbon_input_caspC, mapping = aes(x=open_water_day, y = carbon_mg), color = "#990F26", fill = "#990F26")+
  theme_classic() +
  theme(panel.background = element_rect(fill = "transparent"), # bg of the panel
        plot.background = element_rect(fill = "transparent", color = NA),
        axis.text.x = element_text(size=28),
        axis.text.y = element_text(size = 28), 
        legend.text = element_text(size = 22),
        legend.text.align = 0,
        axis.title = element_text(size = 30), 
        legend.title = element_text(size =22),
        legend.background = element_rect(fill=alpha('white', 0.4)),
        plot.title = element_text(size = 30, hjust = 0.5))+
  scale_x_continuous(breaks = c(-20, 0, 20, 40, 60, 80, 100))+
  scale_y_continuous(breaks = c(0,0.25,0.5,0.75,1,1.25))